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We characterize vibrational motion occurring at low temperatures in dense suspensions ol soft 
repulsive spheres over a broad range of volume fractions encompassing the jamming transition at 
(T = 0, ip = ipj). We find that characteristic time and length scales of thermal vibrations obey 
critical scaling in the vicinity of the jamming transition. We show in particular that the amplitude 
and the time scale of dynamic fluctuations diverge symmetrically on both sides of the transition, 
and directly reveal a diverging correlation length. The critical region near tpj is divided in three 
different regimes separated by a characteristic temperature scale T* (</?) that vanishes quadratically 
with the distance to ipj. While two of them, (T < T*(tp), ip > tpj) and (T < T*(ip), <p < <pj), are 
described by harmonic theories developed in the zero temperature limit, the third one for T > T*(ip) 
is inherently anharmonic and displays new critical properties. We find that the quadratic scaling of 
T*(ip) is due to nonperturbative anharmonic contributions, its amplitude being orders of magnitude 
smaller than the perturbative prediction based on the expansion to quartic order in the interactions. 
Our results show that thermal vibrations in colloidal assemblies directly reveal the critical nature of 
the jamming transition. The critical region, however, is very narrow and has not yet been attained 
experimentally, even in recent specifically-dedicated experiments. 

PACS numbers: 05.10.-a, 05.20.Jj, 64.70.qj 



I. INTRODUCTION 

Important theoretical efforts have recently been de- 
voted to the study of the jamming transition and the 
critical properties of the so-called J-point [l[. An im- 
portant feature of a system approaching the jamming 
transition is the emergence of a peculiar density of states 
which differs from the ones of ordinary solid materials. 
The studies of athermal packings of soft repulsive parti- 
cles above the jamming transition 043 > an< ^ °f nar d par- 
ticles below jamming [8[ both unveiled that the density of 
states acquires an excess of modes at low frequencies, and 
thus increasingly differs from the usual Dcbyc behavior as 
jamming is approached. The presence of an excess of har- 
monic modes which are distinct from usual plane waves 
has been proposed as an explanation to many unusual 
low-temperature properties of disordered materials [l[. 

Motivated by these theoretical results, several recent 
experiments aimed at measuring the vibrational motion 
of a variety of dense colloidal assemblies [i^-flol] . Exper- 
imentally, the normal modes are accessed by diagonaliz- 
ing the displacement correlation matrix to obtain quan- 
titative information about time scales (eigenvalues) and 
length scales (eigenvectors) of the vibrational motion. 
These experimental studies were made possible through 
the development of increasingly accurate tools to record 
single particle motion in colloidal assemblies. Similar 
studies for standard molecular glasses are not possible 
since particle motion cannot be visualized directly in that 
case, but vibrational motion can also be accessed in sys- 
tems such as driven granular media, where attempts to 
characterize modes were recently reported [l6| . 

We wish to characterize the relation, if any, between 
these two sets of efforts. Arc the normal modes deter- 



mined in experimental work related to the ones studied 
theoretically and numerically near the jamming transi- 
tion? What are the necessary conditions to observe ex- 
perimentally the 'low-frequency' or 'soft' modes discussed 
theoretically? What are the real space signatures in sin- 
gle particle trajectories of the anomalous density of states 
that develops near jamming? 

More fundamentally, these questions raise the issue of 
the validity of a harmonic description of the dynamics of 
amorphous materials at low temperature. Indeed, the- 
oretical studies are performed in the zero-temperature 
limit where a harmonic description (or an effective one 
for hard spheres H) is justified. On the other hand exper- 
iments are inevitably performed at finite temperatures in 
colloidal systems, or under gentle vibrations for granular 
materials. Therefore, understanding to what extent fi- 
nite temperatures affect the theoretical predictions and 
quantifying the role of anharmonicity are two important 
problems which are tackled by the present work. 

The fact that finite temperature can be an extremely 
singular perturbation close to jamming can easily be 
grasped recalling that an excess of low energy modes may 
lead to large fluctuations which can eventually destabi- 
lize the system, as it happens for two-dimensional crys- 
tals (l7j . For the jamming problem, there exist two pos- 
sibilities: (1) The harmonic description remains valid at 
low enough temperature. Then, the flat density of states 
would imply arbitrarily large fluctuations and would sug- 
gest the melting of the disordered solid, which would then 
invalidate a purely harmonic description. (2) The har- 
monic description breaks down approaching the jamming 
transition at any value of the temperature. This would 
allow the existence of a disordered solid but would also 
imply that anharmonicity is crucial to prevent the solid 
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from breaking apart. In both cases, the conclusion is 
that the harmonic description must break down at any 
finite temperature close enough to tpj. Because exper- 
imental |a M, El El El and numerical [H-|2l| works 
show that close to jamming a disordered solid is stable, 
the second possibility is clearly the correct one, and an- 
harmonicity must play an increasing role approaching the 
jammin g tr ansition. This issue was recently discussed in 
Rcfs. [5, 22, 23], but our conclusions will differ somewhat 
from these earlier studies. 

In this work we provide a complete description of the 
role of finite temperatures and anharmonicity on the crit- 
ical behavior close to the jamming transition and the 
understanding of the connection, or the lack thereof, be- 
tween vibrational properties studied in experiments and 
theory. We explain when it is correct to assume that ther- 
mal vibrations and eigenmodes of the dynamical matrix 
coincide, and the range of parameters where single parti- 
cle motion becomes dominated by genuine 'soft' or 'low- 
frequency' modes of the type discussed in the context of 
the jamming transition. This issue is also important the- 
oretically, as it could help us to understand better the 
validity of analytical replica calculations where thermal 
fluctuations are taken into account [H HH, but vibra- 
tional eigenmodes are not directly described. 

Finally, because we focus on thermal fluctuations, we 
cannot, in principle, address the more complicated sit- 
uation where particle motion does not result from an 
equilibrium thermal bath, as is the case in granular me- 
dia [IB, 23]. However, the type of 'vibrational heterogene- 
ity' that we reveal through the study of thermal vibra- 
tions is strongly reminiscent of the findings reported in 
two distinct experimental studies of dynamic heterogene- 
ity in driven granular media across the jamming tran- 
sition (26l . [27j . In particular, we obtain at finite tem- 
perature a nonmonotonic dependence of both a dynamic 
susceptibility and a correlation length scale very simi- 
lar to the experimental findings for dissipative grains. 
They stem, in our case, from the existence of a corre- 
lation length that is divergent when the temperature is 
strictly zero. Surprisingly, this divergent length scale is 
physically distinct from the 'isostatic' length scale which 
is central to the description of the density of states near 
jam ming or the response function of athermal pack- 
ings [231. We will connect this last finding with recent 
investigations of random networks Q. 

This paper is organized as follows. In Sec. [TT| we de- 
scribe our numerical model and methods. In Sec. IIII1 
we present numerical results for the dynamic behavior of 
low-temperature packings of soft repulsive spheres across 
the jamming density. In Sec. IIV1 we analyze more pre- 
cisely the scaling behavior of these dynamic quantities, 
and the emergence of diverging time scales and length 
scales near jamming. In Sec. [V] we show that these di- 
vergences result from the existence of diverging dynamic 
fluctuations and of a diverging correlation length scale. 
The finite temperature critical behavior is analyzed in 
Sec. I VII In Sec. IVIII we numerically determine the tem- 



perature scale for the emergence of anharmonicity, and 
provide a theoretical analysis of its physical content. Fi- 
nally, we use our findings to discuss experimental results 
in Sec. IVnTl 



II. MODEL AND NUMERICAL METHODS 

We use molecular dynamics simulations to study the 
thermal vibrations of sphere packings near the jam- 
ming transition. We focus on monodisperse harmonic 
spheres [29j , interacting through a pairwise potential, 

v(r ij ) = ^(l-r ij /a) 2 e(a-r ij ), (1) 

where Q(x) is the Heaviside function and a is the par- 
ticle diameter. We use system sizes between N = 120 
and N = 64000 particles imposing periodic boundary 
conditions, changing the volume V to adjust the pack- 
ing fraction tp = ira 3 N/(6V). The main results are re- 
ported for N = 8000, while we use larger and smaller sys- 
tems for specific purposes to be specified below. We use 
Newtonian dynamics, and integrate Newton's equation 
of motion in the microcanonical ensemble after proper 
thermalization of the system at the desir ed tem perature. 
The results are reported in units of ^/ma 2 /e for time 
scales, of a for length scales, and of e/fcs for tempera- 
tures, where ks is the Boltzmann constant. 

To produce a low temperature amorphous packing of 
harmonic spheres in three dimensions at volume fraction 
tp, we first perform an instantaneous quench from T = 
00 to T = 10 ~ 5 (well below the glass temperature (l8j |) 
at a large volume fraction, tp = 0.70. We then let the 
system relax and rapidly settle in a metastable state. 
At these very low temperatures (compared to the glass 
temperature) aging effects rapidly become negligible and 
the system simply performs vibrations near a well-defined 
energy minimum. 

We then explore the (tp, T) phase diagram of these 
packings by decreasing T and tp to the desired values. 
Using this method we can follow the same jammed con- 
figuration with different particle overlap or gap when 
changing temperature or density. For each independent 
quench, the jamming density tpj is therefore well-defined 
and can be measured very accurately |30ll3ll |. for instance 
by decompressing the packing at T = towards tpj and 
measuring the volume fraction where the pressure van- 
ishes. Following previous work [2p| . we mainly discuss 
the vibrational behavior of a single, large packing with 
N = 8000 and for which tpj ~ 0.6466. We provide ad- 
ditional data for different packings where necessary. We 
checked that the results reported in this work do not de- 
pend on the specific choice of a packing. 

Note that all dynamic observables discussed in this 
paper are measured after aging effects have terminated, 
and are thus stationary or independent of the waiting 
time, even though the dynamics is not ergodic. We 
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merely explore the statistical properties of a representa- 
tive metastable state. Thus, our studies are qualitatively 
distinct from the non-stationary aging dynamics studied 
in Ref. @, and from the athermal studies of Refs. (32| — 
[34[ where an externally imposed shear flow affects the 
structure and dynamics of the packings. 

III. DYNAMICAL BEHAVIOR AT FINITE 
TEMPERATURES NEAR JAMMING 

A. Mean-squared displacement 

In order to characterize vibrational motion, we focus 
on single particle dynamics. The most basic correlation 
function is the mean-squared displacement (MSD), which 
is easily measured experimentally in colloidal systems, 

1 N 

A oW = ^E(i A ^wi 2 )' ( 2 ) 

1=1 

where Arj(t) = rj(t) — 1^(0), and rj(t) represents the 
position of particle i at time t. The brackets in Eq. @ 
represent an average over long trajectories performed at 
a given T and <p (but nevertheless restricted to a single 
metastable state). We will be particularly interested in 
the long-time limit of the MSD, which we call the Debye- 
Wallcr (DW) factor. 

Near the jamming transition, there is a slight techni- 
cal problem arising in these measurements because the 
packings produced numerically contain a small fraction 
of 'rattlers', that have peculiar dynamic properties since 
they are less connected than most particles and therefore 
have distinct dynamical properties. An annoying feature 
is that rattlers move much more than the other parti- 
cles, and so they might completely dominate the sum in 
Eq. ((2J, thus yielding results which do not reflect the 
representative behavior of the system. 

We used two distinct methods to prevent rattlers from 
spoiling the measurements. A first, simple method con- 
sists in computing a modified version of the MSD, 

1 N 3 
ALd (*) = jj E (|A ri (i)|-2) ' (3) 

which gives nearly no statistical weight to the particles 
that move much further than the average. The prefac- 
tor 3 in the definition © is such that A^ od (t) = Ag(i) 
when the distribution of single particle displacement is 
Gaussian. In a second method, we directly remove the 
rattlers and measure 

N' 

A 2 (t) = — £(|Ar 4 (t)| 2 ), ( 4 ) 

i=l 

where N' is the total particle number which are not rat- 
tlers. Although rattlers can be defined unambiguously in 



the jammed regime, ip > ipj, by measuring contact num- 
bers in the ground state at T = 0, we cannot apply this 
method to the unjammcd regime, ip < <pj. We have ap- 
plied the following procedure which can be used in both 
regimes. We first measure the DW factor of each particle. 
We then calculate the median value, which is less affected 
by rattlers than the averaged value. Then we regard the 
particles whose DW factor is 5 times larger than the in- 
termediate value as rattlers. We have confirmed that this 
procedure and the analysis at the ground state give the 
same result in the jammed regime. We have also con- 
firmed that the definitions Eqs. ([3]) and (UJ) yield similar 
results at all densities and temperatures. We note that 
definition <j3j> is conceptually extremely simple and does 
not require any empirical criterion to identify rattlers, 
and could therefore be very useful for the analysis of ex- 
perimental data. In the following, we present the data 
obtained from the definition (QJ. 

In Fig. [TJ we present typical results for the behavior of 
the MSD obtained at constant temperature, T = f0~ 8 , 
and a range of volume fraction encompassing the jam- 
ming density ipj. The time dependence of the MSD is 
typically characterized by two distinct time scales. At 
very short times, the MSD is purely ballistic and fol- 
lows A 2 (t) — Tt 2 , which simply results from the short- 
time integration of Newton's equation of motion for par- 
ticles thermalized at temperature T. Deviations from 
this trivial behavior occur after a microscopic time scale, 
To, which corresponds to the time where particles start 
to feel the influence of their neighbors. This is indicated 
by the open squares in Fig. [TJ We shall characterize tq 
more accurately in Sec. lIV Al below. but the data in Fig.Q] 
indicate that To decreases with <p. 

For times t > tq , the time evolution of the MSD crosses 
over to an intermediate time regime where it obeys dif- 
fusive behavior, before saturating at long times to a 
plateau, the DW factor A 2 (oo), which decreases rapidly 
upon compression. The approach to the plateau value 
occurs after a typical time scale, t*, marked by the filled 
squares in Fig. [TJ The precise definition of t* is also dis- 
cussed below. The data indicate again that t* decreases 
strongly when the system is compressed. The fact that 
the DW factors remain finite in the long-time limit shows 
that our systems are genuine amorphous solids [35| . and 
that they do not melt despite having a flat density of 
states at T = and ip = ipj. 

We now extract the Debye- Waller factor from the long- 
time limit of the mean-squared displacements. We gather 
our results for a broad range of densities and tempera- 
tures in Fig. [TJ We find that the DW factor decreases 
upon compression at fixed temperature, mirroring the 
behavior of both time scales To and t*. When the sys- 
tem gets denser, the particles have less space to move, 
their vibrations are more constrained, and so the typical 
amplitude of the vibrations and the related time scales 
become smaller. Therefore, a smooth evolution of dy- 
namic quantities with density can be anticipated from 
trivial reasons, unrelated to the specific physics of the 
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FIG. 1: Top: Time dependence of the mean-squared displace- 
ments (MSD) at constant temperature, T = 10 -8 , and var- 
ious volume fractions increasing from top to bottom across 
the jamming transition. Open squares indicate the micro- 
scopic time scale to, filled squares indicate the long timescale 
t*, which marks roughly the convergence of the MSD to its 
long-time limit. Both time scales decrease with tp, but their 
ratio is maximum near tpj. Bottom: Volume fraction depen- 
dence of the Debye- Waller (DW) factor (long-time limit of 
the MSD) for different temperatures. The T — jamming 
singularity at tpj ~ 0.6466 (vertical line) strongly influences 
the DW factor when T becomes very small, T < 10 -7 , and tp 
is very close to the critical density, 0.64 < tpj < 0.655. 



jamming transition. 

Two interesting features emerge from the DW data 
in Fig. Q] Firstly, the scaling with temperature is dif- 
ferent on both sides of the jamming transition. This is 
not surprising, as many other observables have a similar 
dual behavior (2f|, which corresponds to the two differ- 
ent situations obtained when T — >• 0. While the system 
corresponds to a hard sphere assembly when if < tpj, it 
is instead a jammed soft solid for tp > tpj. Accordingly, 
the DW factor converges to its hard sphere value below 
jamming and remains finite as T — > 0, while it vanishes 
linearly with T above jamming, as in ordinary solids pTj . 

Secondly, we also note that as T becomes smaller, the 
volume fraction dependence of the DW factor becomes 



more and more singular in the vicinity of tpj, which sug- 
gests that some singularity will emerge in the T — > 
limit, to be discussed below. Here, we simply empha- 
size that for the harmonic sphere system, temperature 
must be at least lower than T = 10~ 7 (in units of the 
spring constant e) to notice the emergence of a singular 
volume fraction dependence of the DW factors near the 
jamming transition. For larger temperatures, the DW 
factor smoothly depends on T and tp as a result of the 
compression but the jamming singularity does not man- 
ifest itself. Anticipating the discussion of experimental 
work to come in Sec. IVIIIi we also remark that a singu- 
lar behavior can be seen in the data when the DW factors 
take extremely small values, typically A 2 (oo) < 10 -4 (in 
units of the particle size) , which corresponds to resolving 
distances with a precision greater than c/100, which can 
be technically quite challenging since it corresponds to 
about lOnm for colloids with a = 1/xm. 



B. Velocity autocorrelation and density of states 

Another way of looking at single particle motion is to 
record the velocity autocorrelation function defined by 



d(t) 



1 



N' 



3N'T 



(5) 



where Vj(i) denotes the velocity of particle i at time t. 
Note that we again remove rattlers from the definition in 
Eq. ([5]). When Fourier transformed, this becomes what 
we (abusively) call the 'density of states' (DOS) [3t| : 



d(w) 



dt cos(tdt)d(t). 



(6) 



Note that d(u>) can be measured arbitrarily far from the 
harmonic limit and so it does not represent a genuine den- 
sity of states, but it does reduce to the true DOS (as mea- 
sured by diagonalizing the dynamical matrix) when the 
harmonic approximation holds [23j | . Below, we will care- 
fully determine when this assumption is correct. Note 
also that by definition one has 
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which shows that both types of measurements (density of 
states and mean-squared displacements) are directly re- 
lated and, in fact, contain nearly equivalent information. 
This remark suggests that, when carefully analyzed, the 
MSD themselves should directly reveal the anomalous 
features of the DOS, as we shall demonstrate. 

In Fig. [21 we present the frequency dependence of d(co) 
at the same set of temperatures and densities as in Fig. [1] 
At low frequency, one gets the expected Debye behavior, 
d{oS) ~ oj 2 for our three-dimensional systems. Moving to 
larger frequencies, we see that the DOS exhibits a single 
peak below jamming, tp < tpj, and this peak broadens 
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and [2] confirms that both quantities contain equivalent 
information in the sense that they both reveal that vibra- 
tional dynamics is controlled by two distinct time scales, 
To = tt/ujq and t* = Tt/u*. The former indicates when 
particles start to feel their neighbors and deviations from 
ballistic motion appear. The latter corresponds to the 
convergence of the MSD to its long-time limit, and also 
to a crossover to a Dcbye behavior in the DOS, meaning 
that interesting dynamics arises in the range To t <C t* , 
or to* -c oj -C wo in the frequency domain. 

Note also that the DOS reported in Fig.[3]arc fully con- 
sistent with previous determination in both jammed Q 
and unjammed regimes Q. The main feature noticed in 
previous work in these DOS data is the emergence of the 
low-frequency crossover u*, which we can alternatively 
interpret, in the time domain, as the time it takes for the 
MSD to reach its plateau value, see Fig. [TJ 

We extract t* = ir/ui* for a broad range of volume 
fractions and temperatures, see Fig. [5J As in the case of 
the DW factor, t* shows very different temperature de- 
pendence in jammed and unjammed regimes. It is nearly 
independent of temperature in the jammed regime but 
scales as t* ~ y/T, reflecting the fact that temperature 
trivially renormalizes the unit time scale in hard sphere 
systems. As temperature gets lower, we observe that t* 
acquires a strong density dependence very close to ipj, 
changing by several orders of magnitude over a narrow 
density range. By contrast, t* changes smoothly when T 
is not low enough, as a consequence of the compression. 



FIG. 2: Top: Frequency dependence of the density of states 
at constant temperature, T = 1CP 8 , and densities increas- 
ing from top to bottom. Parameters are as in Fig. [1] Open 
squares indicate the microscopic frequency tJo = tt/to, while 
filled squares indicate the low-frequency crossover oj* = n/t* . 
Bottom: Volume fraction dependence of t* = n/uj* for differ- 
ent temperatures. The vertical line indicates ipj. 



with increasing density and transforms into a plateau in 
the jammed regime, <p > tpj. To determine the position 
of this low-frequency characteristic behavior over the en- 
tire density range, we define the crossover frequency w* 
as the peak position of the quantity d(uS)/u>. In Fig. [5] 
we indicate the frequency uj* with filled squares, which 
appear as good indicators of the low-frequency crossover 
towards the Debye power law. Clearly, lu* shifts to larger 
frequency when (p increases. 

At large frequency, d(tu) undergoes a second crossover 
towards either an algebraic d(uS) ~ a;" 2 behavior below 
jamming, or towards an even steeper decay above jam- 
ming. This behavior simply mirrors a qualitative change 
in the short-time behavior of the velocity autocorrela- 
tion function from exponential to Gaussian over the same 
density range (37|. We report the microscopic frequency 
luq = tt/tq as open squares in Fig. [2] which confirms that 
To governs the high-frequency crossover in the DOS. 

The comparison between the MSD and DOS in Figs.Q] 



IV. CRITICAL BEHAVIOR 

In the previous section, we discussed the qualitative 
behavior of various time scales characterizing the vibra- 
tional motion of dense packings of harmonic spheres. In 
particular, we noticed that singular density dependences 
seem to emerge when temperature decreases. In this sec- 
tion, we focus on these low-temperature singularities and 
show that they correspond to zero-temperature algebraic 
divergences rounded by finite temperature effects, and 
therefore to an underlying critical dynamics. 



A. Elementary time scales and length scales 

From the data reported in Figs. Q] and [2j one finds that 
the dynamics shows an abrupt crossover at low temper- 
ature when ipj is crossed. Without further analysis, it is 
however not obvious to decide when a nontrivial critical 
behavior becomes manifest, as compression itself influ- 
ences the dynamics. This suggests that one should first 
rcnormalize time scales and length scales in terms of the 
microscopic ones. In the case of the jamming transition, 
compared to usual critical phenomena, this is particularly 
important since the dynamics below and above jamming 
are very different: dynamical behaviors are trivial for 
ip <C j and if 3> tpj but in a quite distinct way. 
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FIG. 3: Volume fraction dependence of the microscopic time 
scale to, Eq. (JSJ), at various temperatures decreasing from 
bottom to top. The vertical line indicates ipj. These data 
were used to draw the open squares in Figs. [T] and [21 Inset: to 
is replotted as a function of the scaling variable \tp — tpj\/yT. 
Upper and lower dashed lines correspond to T — > limits for 
hard and soft spheres, respectively, Eqs. (J9] ITOjl . 

Our first step is to analyze the short-time dynamics 
in order to understand the behavior of the microscopic 
time scale tq discussed above, which also sets the high- 
frequency cutoff in the DOS d(u>). In practice, we mea- 
sure To by analyzing the short-time decay of the velocity 
autocorrelation function d(t) in Eq. ([5]), and quantita- 
tively determine r from the definition 

d(t = r ) = 1/e. (8) 

The data for r are shown in Fig. [3] These data were also 
used to determine the position of the open squares in the 
MSD data of Fig. Q] and in the DOS data of Fig. H 

The short-time behavior of the velocity autocorrelation 
function d(t) is very different in jammed and unjammed 
regimes at low temperature. In the zero-temperature 
limit, d(t) display an exponential decay when tp < tpj 
due to two-body collisions in the hard sphere limit, whose 
characteristic time scale can be deduced from Enskog the- 
ory (roughly the time it takes to a particle to collide bal- 
listically with a neighbor) [33] , 

1 ^ 8pg{*+) p VT 
t 3 V m \(p-(pj\' 

where the final approximation only holds very close to 
jamming. We stress that this two-body collision time 
scale decreases very rapidly as tpj is approached from 
below because the typical gap between neighboring par- 
ticles also vanishes rapidly, scaling as (tp — tpj). This is 
quantitatively taken into account in Eq. © through the 
contact value of the pair correlation function, g(o~ + ). 

In the low temperature limit above jamming, tp > ipj, 
d(t) shows instead a Gaussian decay, whose characteris- 
tic timcscale can now be estimated as the inverse of the 



Einstein frequency [l7l ]. 



— ~ J— / drg(r)\7 2 v(r) ~ const, (10) 
t y 3m J 

where g(r) is the pair correlation function, and the final 
approximation means that To has no singular dependence 
upon either T or tp near jamming. 

The data in Fig. [3] indicate that To interpolates 
smoothly between the two limits discussed above on both 
sides of the transition at finite temperatures. The main 
panel confirms that it is proportional to VT below jam- 
ming, and becomes independent of temperature in the 
jammed regime. As in the case of t*, it is only in the 
very low temperature limit that a singular density de- 
pendence emerges. In the inset of Fig. [3l t is replotted 
as a function of the rescaled variable x = \tp — pj\/VT, 
in order to reconcile the behavior on both sides, as is 
the case for most observables near jamming (Til. l25l. l32| . 
Indeed we get a very good collapse of all the data along 
two distinct branches, confirming the simple dependences 
discussed above. We emphasize that even though t is 
controlled by relatively simple physics and it simply sets 
the microscopic time scale for the dynamics, it changes 
rather dramatically near jamming at low temperatures. 

Having identified the elementary time scale To on both 
sides of the jamming transition, we shall now focus on 
the elementary length scale. Below jamming, as dis- 
cussed above, particles move ballistically over the typical 
interparticle distance in a time To- Thus, the elementary 
length scale obtained by ballistic motion reads: 

i = VTt . (11) 

The same relation holds on the other side of the jam- 
ming transition. In that case the elementary length 
scale is given by the typical displacement due to a har- 
monic mode with Einstein frequency (|10[) . leading again 
to Eq. (fTTj). In both cases, £q represents the distance 
travelled ballistically over the microscopic time scale t . 

In conclusion, in this section we have identified the el- 
ementary length scales and time scales, Eqs. (|9l [TOl [TT1) . 
which characterize noncollective dynamics. In the next 
section we shall show that approaching the jamming tran- 
sition, there emerge time scales and length scales which 
actually become much larger than To and £q, thus sig- 
nalling the existence of bona fide collective dynamics. 

B. Critical behavior of adimensional data 

We now reconsider the dynamical behavior studied in 
Sec. IIIII and analyze the data obtained by renormalizing 
all observables in terms of the elementary time scales and 
length scales. 

We define a rescaled Debye- Waller factor as 

A oc (tp 1 T) = — p — = . (12) 
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FIG. 4: Critical behavior for rescaled dynamic observables. 
Top: Rescaled Debye- Waller factor, Eq. H2|. Bottom: 
Rescaled time scale for vibrational motion, Eq. (|13[) . In 
both panels, dashed lines represent the critical divergence 
~ \<p — ipj]^ 1 ^ 2 , Eq. (|18|) . and the vertical line indicates tpj. 



and we similarly renormalize the long time scale t* in 
terms of tq : 



TO 



(13) 



A visual interpretation of these adimcnsional quantities 
can be given from the symbols in Fig. [TJ since Aoo and 
t respectively quantify the vertical and horizontal dis- 
tances between the filled and open squares, i.e. between 
the short-time and long-time dynamics. 

The numerical data for these two quantities are re- 
ported in Fig. SJ Clearly, they both arc smooth func- 
tions of the density when temperature is not very low 
(but still much smaller than the glass temperature), e.g. 
T = 10~ 5 , and T = 10 ~ 6 , and they also are smooth func- 
tions of the temperature when tp is not in the immediate 
vicinity of tpj. This implies that in this regime, the 'low- 
frequency' crossover lo* is not very 'low', i.e. not strongly 
decoupled from the natural microscopic time scale To- In 
the similar vein, in the non-critical regime the DW factor 
is not very much larger than its natural scale, and thus 



thermal vibrations are neither 'anomalous' nor charac- 
terized by 'soft' modes which would reflect an underlying 
nontrivial or collective dynamics. 

The situation becomes more interesting at very low 
temperature, T < 10~ 7 and for densities near jam- 
ming, 0.64 < ip < 0.655, because adimcnsional quantities 
may become very large, and additionally they acquire 
a nonmonotonic density dependence across <p j, strongly 
suggesting the existence of divergences as T — > on 
both sides of the transition. Thus, we conclude that 
the anomalous 'low-frequency' behavior of the DOS ap- 
proaching jamming corresponds to the appearance of an 
adimcnsional timescale r that becomes large, t» 1, and 
an 'anomalously' large adimcnsional amplitude of the vi- 
brational motions, Aoo 3> 1. Thus, when analyzed from 
the viewpoint of thermal vibrations, the jamming transi- 
tion is accompanied by a critical slowing down which wc 
now confirm through a scaling analysis. In Sec. [Vj wc 
will show that this slowing down is also accompanied by 
a diverging correlation length. 



C. Critical behavior and scaling in the T — > limit 

We now connect the two critical behaviors found above 
and estimate the corresponding T — > divergences using 
a scaling analysis. This analysis amounts to assuming 
that slow vibrational motion is controlled by a single 
dominant timescale, t* , i.e. that the DOS at low fre- 
quency has a non-trivial scaling behavior of the form 



d{u) = T f± f-^J , w < 1/t , 



(14) 



where f± (x) corresponds to the universal behavior of the 
density of state below (— ) and above (+) jamming, re- 
spectively. We checked that numerical data are consis- 
tent with this scaling assumption. Both functions f±{x) 
are proportional to x 2 for x — > in order to recover the 
usual Debye law. For x — > oo, f+(x) tends to a con- 
stant to recover the flat part of the density of states, 
while / (x) decreases to reproduce the hard sphere be- 
havior [50(. 

By definition of the MSD and of the DOS we get 

A 2 (<) = 6T r dw^[l - cosfwt)]. (15) 
Jo ^ 

In the long-time limit, the scaling behavior in Eq. (fT^Ij) 

yields 



rp />OC 

A 2 (oo) w ±H / dx 
which simply implies that 



/±(*) 



A- 



(16) 



(17) 



This result unveils that, to leading order, the data in both 
panels of Fig. U are in fact equivalent, up to a numeri- 
cal prcfactor stemming from the difference in the scaling 
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functions f±(x) expressed by Eq. (|T6]) . The result in 
Eq. ([T7| has a simple physical interpretation, because it 
says that a large amplitude for vibrational motions devel- 
ops near jamming, and that the time it takes for particles 
to fully explore these vibrational degrees of freedom also 
becomes larger, the scaling between (adimensional) times 
and lengths obeying a simple diffusive scaling. 

Finally, by using literature results concerning the scal- 
ing of uj* approaching the jamming transition at T = 
jH, Q , we find that r diverges symmetrically on both 
sides oi <pj &s \ip — (pj\~ x / 2 . For our rcscalcd variables, 
this implies that when T — > 

a^-t- i^-^r 1 / 2 , (is) 

on both sides of the jamming transition. These scaling 
laws are in good agreement with our data, as shown with 
dashed lines in Fig. |4j We have not tried to measure 
the exponent in Eq. (fT5)) directly to detect a possible 
deviation from the value —4, but our numerical results 
suggest that this correction, if it exists, is quite small, in 
agreement with recent numerical work [§], [5l[ . 

As is clear Fig. 21 the divergence is cut off by using a fi- 
nite temperature. As we shall discuss later, this is due to 
the appearance of anharmonicity at finite temperatures. 

V. DIVERGING CORRELATION LENGTH 

In the previous section, we found that the amplitude 
of the vibrational motion, and the associated time scale, 
may become nontrivial, i.e. much larger than the ex- 
pected elementary values corresponding to non-collective 
dynamics. It is therefore natural to expect some kind of 
diverging correlation length scale. It is the purpose of 
this section to unveil and study this correlation length. 

A. Four-point correlation function 

Since we identified the mean-squared displacement as 
the quantity of interest to pinpoint the existence of the 
transition, we now seek spatial correlations between the 
particle displacements during the vibrational motion. 
This amounts to measuring the extent of spatial correla- 
tions in the dynamics of glassy states. We can therefore 
use the machinery developed to study dynamic hetero- 
geneity in glassy materials jJ{-> . First, we measure the 
following four-point susceptibility 

X4 (t)=iV[( M (i) 2 )-(^)) 2 ], (19) 

defined as the variance of the spontaneous fluctuations of 
the instantaneous value of the volume-averaged mobility 
/i(i), defined as 



Note that the thermal average in the denominator of 
Eq. (f2T)|) does not include an average over particles, and 
thus particles with larger displacements do not dominate 
the sum in the definition of n{t). This implies in par- 
ticular that rattlers are not dramatically affecting the 
numerical measurements of Xi(t)- 

While Xi(t) measures the total amplitude of spatial 
correlations, direct access to a correlation length scale is 
provided by the spatially resolved correlator, which we 
study in the Fourier space: 

1 N 

S4(q,t) = jj ^(|^(tK(t) e ^|), (21) 

3,k=l 

where rjfe = r_, (0)— rfc(O) is the separation between parti- 
cles j and k at t = and [ij (t) is the contribution of par- 
ticle j to the mobility, fj,j(t) = \Ar j(t)\ 2 /{\Ar jt)\ 2 ) - 1. 

Two remarks before showing our results. First, we have 
included all particles in the analysis of the dynamic het- 
erogeneity, because the rattlers do not dominate the sum 
in Eqs. (fT9l l2Tj) . contrary to the MSD. We have numeri- 
cally checked that our results in this section are identical 
if we remove the rattlers. Second, our definition of X^iP) 
does not explicitly involve a 'probing' length scale as in 
other works [111 [38|. In fact, using the mobility defined 
in Eq. ([2U| is equivalent to selecting automatically the 
appropriate probing length at each time (26|. 

We find that Xi{t) increases with t at short times, 
but it does not display a clear maximum at intermediate 
times as is usually found in systems near a glass tran- 
sition [39], but saturates to a long-time limit discussed 
below. Here we wish to characterize the spatially hetero- 
geneous dynamics relevant for Aqo and r. We therefore 
concentrate on time scales of the order of t* discussed 
above, and we measure 

Xi&,T) = X4(t = t*), (22) 

and the corresponding structure factor: 

S i (q;<p,T) = S i (q,t = t*). (23) 

Representative numerical results are shown in Fig. [5] 
As for the adimensional quantities Aqo and r in Fig. |4j 
X4 also shows a striking nonmonotonic density depen- 
dence across ipj, which becomes more pronounced as T 
decreases, and would seem to diverge on both sides of the 
jamming transition as T — > 0. These results strongly sug- 
gest that the slowing down and excess amplitude of ther- 
mal vibrations are the result of collective, spatially corre- 
lated dynamics. Remark again that Xi acquires a strong 
nonmonotonic density dependence in the close vicinity of 
95,7, but only for temperatures which arc at least below 
T ~ 10~ 7 , while it remains nearly constant and rather 
insensitive to the underlying jamming transition every- 
where else. Finally, remark that the critical behaviour 
of xa emerges in Fig. [5] without further analysis or nor- 
malization by an elementary 'correlation unit' (5l| . This 
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FIG. 5: Top: Volume fraction dependence of the dynamic 
susceptibility x*(t — Eq. (|22l) . for various temperatures 
showing the emergence of 'vibrational heterogeneity' at low 
enough T close to tp j (shown with the vertical line) . Bottom: 
Four-point correlator S4(q,t = t*) at T = 1CF 8 and various 
volume fractions exhibits a nonmonotonic density dependence 
(note the nonmonotonic evolution of the ip labels) . 



makes it a very relevant quantity to be measured exper- 
imentally, as it directly reveals the dynamic criticality 
associated to the jamming transition. 

We now focus on spatial correlations to extract the 
corresponding correlation lengh scale over which thermal 
vibrations occurring over a time scale t* are correlated. 
In Fig.O we present the evolution of 54(g) for T = 10~ 8 
across the jamming density. As suggested by the behav- 
ior of X4, we find that 54(g) is also characterized by a 
nonmonotonic behavior with ip, spatial correlations be- 
ing modest at both low and large 93, and being maximum 
near ipj. Note that 54(g — > 0) represents the volume 
integral of the spatial correlation between particle dis- 
placement, and it is thus natural that the low-wavevector 
limit of 54(g) follows the same trend as X4 HO]- However, 
the spatially resolved correlator allows us to directly ex- 
tract a correlation length through the q dependence of 
54(g). Clearly, the approach to the low-g plateau occurs 
for lower g when ip w ipj, which indicates that the cor- 



responding correlation length is also maximum there at 
finite T, and would diverge at fj in the T — » limit. 

Before discussing the divergence at T = in more de- 
tail, we wish to emphasize the striking similarity between 
the numerical results displayed in Fig. [5] quantifying the 
the 'vibrational heterogeneity' of our thermally excited 
packings of soft particles, and a similar nonmonotonic dy- 
namic heterogeneity reported for gently vibrated grains 
in Rcfs. 



27j . Although the excitation mechanisms 



and the physical systems are quite distinct, it is very 
tempting to suggest that both sets of measurements are 
actually probing the same underlying dynamic criticality 
associated to the jamming transition. 



B. Critical behavior and scaling in the zero 
temperature limit 

We now connect these observations of spatially hetero- 
geneous vibrational dynamics with the critical behavior 
of t and Aqo, extending the previous scaling analysis. 
Our starting point to evaluate \4 is to assume that fluc- 
tuations are Gaussian and, hence, four-point correlation 
functions can be written as products of two-point func- 
tions. This is correct within the harmonic approxima- 
tion, which is valid in the zero-temperature limit we are 
considering for this analysis. With this assumption, we 
can obtain an analytical expression of the long-time limit, 
t — > 00 , of the four-point susceptibility as 



T 



A 2 (oo) 



du> 



d{uj) 



(24) 



Here w m i n is the lowest frequency of the density of state 
for a given system. When the system size becomes large, 
w m i n ~ L^ 1 — > and d(cu) cx oj 2 for small oj, and then, 
the integral displays an infrared divergence which is cut- 
off at a system size dependent value, as we have numer- 
ically confirmed. This is because the lowest frequency 
mode is a plane wave for which the motion of all the 
particles in the system are fully correlated [39| . 
We focus instead on finite times, Xi ~ Xi(t = 



x*(f) 

from which we deduce 
X4 ^ 



T 



A 2 (i*) 



d{uj) 



d(ui*)u)* 



(25) 



(26) 



Finally, if we assume (to be confirmed in a short while) 
an Ornstein-Zcrnike form for the correlation length we 
also obtain the relation between susceptibility and length 
scale, £4 - y/xi. 

To summarize our scaling analysis, we have shown that 
the following quantities are all connected in a simple way 
at the level of their scaling behaviors, and we obtained 
the following critical scalings: 



X4 ~ £4- 



(27) 
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FIG. 6: Test of the scaling behavior of S4,(q) using the as- 
sumption that £4 ~ x/Aqc. Note that this data collapse in- 
volves no free parameter since S4 and Aoo result from inde- 
pent measurements. Dashed lines are fits to the Ornstein- 
Zernike functional form, Eq. (|28ll . 



Note that in order to derive these results we only used 
the harmonic approximation (valid when T — > 0) and 
the scaling assumption (j!4j) . If we add the information, 
obtained numerically in Refs. [3, HI, that r diverges as 
\f — (pj]^ 1 / 2 then we also get the zero temperature di- 
vergence of Aqo oc \(p — tpj]^ 1 / 2 and of £4 cx \<p — ipj]' 1 / 4 
approaching the jamming transition. As discussed before 
and is clear from our data, any finite temperature cuts 
off these divergences. 

According to our analysis, the diverging correlation 
length is proportional to the square root of the diverging 
excess amplitude of the thermal vibrations, £4 ~ V^-oo- 
We directly test this hypothesis in Fig. [5] where we 
present S4(q)/A ao as a function of q^/A^. Clearly, all 
the data are collapsed very well along two branches cor- 
responding to both sides of the jamming transition. We 
emphasize that this data collapse involves no adjustable 
parameter, but makes use of two independent sets of mea- 
surements. Additionally, we show that both branches are 
correctly described by an Ornstein-Zernike wavevector 
dependence, namely 



Si(q) 



c±A c 



1 + c ± q 2 A c 



(28) 



with c± a numerical prefactor which is distinct for the 
two branches. These data directly confirm the validity of 
the analysis in Eq. (|27p. 

We conclude that the appearance of 'anomalously' slow 
and large vibration motion is also associated with spa- 
tially heterogeneous dynamics, which is a form of dy- 
namic criticality |4lj ]. These observations indicate that 
the jamming transition appears as a genuine critical tran- 
sition with both a critical slowing down and a diverging 
length scale. Since the slowing down occurs for thermal 
vibrations, this suggests that a convenient order param- 
eter to describe the transition is the vibrational motion 



within cages deep into the glassy phase. It is intrigu- 
ing that all the critical exponents describing time scales, 
amplitude of vibration, and dynamic heterogeneity are 
related in a trivial manner, Eq. (f2~T)) . Interestingly these 
divergences are symmetric and hold quantitatively for 
both hard and soft spheres, provided adimensional quan- 
tities are considered. 

We note that the length scale quantifying the spatial 
correlations in the excess amplitude of the vibrations di- 
verges, for the harmonic potential used in this study, as 



-1/4 



(29) 



A length scale diverging with the same exponent was al- 
ready discussed in different contexts. In Ref. [3], a diverg- 
ing length is identified from the transverse structure fac- 
tor of the normal modes obtained for soft spheres above 
jamming. This length scale also appears in the analysis 
of Ref. [6| of the heat transport in jammed solids. Finally, 
the same length scale also shows up in a mean-field treat- 
ment of the localization length of the elastic response of 
model networks Q- To our knowledge, for hard spheres, 
this length scale was not discussed before. Surprisingly, 
in both soft and hard cases, the isostatic length scale in- 
troduced in Ref. @ does not directly appear in any of 
the physical quantities we discuss in this work. Finally, 
£4 in Eq. (|2T)|) seems to differ also from the correlation 
lengths measured in Refs. HI El 53]. 

Remark that in our scaling analysis, arguments in- 
volving the average contact number are not directly em- 
ployed. It is only to obtain the divergences of r, Aoo, and 
£4 with \ip— ipj\ that we made use of the scaling of oj* with 
\ip — ipj\ which was related to the distance to isostatic- 
ity @- We stress that none of those divergences is fully 
understood theoretically from a microscopic viewpoint. 
In fact, isostatic arguments [3] trace back the behavior of 
oj* with \<p — tpj\ to the one of the contact number, which 
is not directly derived, but measured from numerical sim- 
ulations. An alternative theoretical approach, which is 
fully microscopic, uses statistical mechanics and replica 
calculations to study the jamming transition (24|. When 
applied to harmonic spheres [13, [25[ , it correctly predicts 
that the jamming transition is characterized by a diverg- 
ing amplitude of thermal vibrations: A M ~ \<p — tpj\~ . 
However, this predicted divergence is only qualitatively 
correct, and the predicted critical exponent is wrong, see 
Eq. (|27|) . It would be interesting to understand better the 
source of this discrepancy from the replica viewpoint 44 1 . 
Our results also suggest that the theory should be recon- 
sidered in order to make precise predictions for a diverg- 
ing correlation length and dynamic susceptibility associ- 
ated to the jamming transition. 



VI. FINITE TEMPERATURE CRITICAL 
BEHAVIOR 

In the following we will fully describe the critical be- 
havior in the (</?, T) phase diagram at the jamming transi- 
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tion. Wc will first analyze the low-T harmonic regimes in 
both jammed and unjammed phases. We shall show that 
their domain of validity is restricted to a temperature 
regime that shrinks approaching ipj. This corresponds 
to an emergent temperature scale, T*(ip), that vanishes 
quadratically with \ip — ipj\. Finally, we shall focus on 
the finite temperature critical regime, corresponding to 
T* « T « 1, that is related to the critical divergences 
when T — > for ip = tpj. 



A. Regime I: Harmonic description above jamming 
and its domain of validity 

In the following we shall establish that approaching 
the jamming transition from the jammed phase, the har- 
monic description used in the previous section is valid 
and allows one to explain the behavior of physical ob- 
servables related to vibrational motions if temperature 
is low enough. Moreover, we shall find the tempera- 
ture scale, T*(ip), above which the harmonic description 
breaks down and the low temperature physical behavior 
changes. 

To this end, we first focus on the temperature depen- 
dence of the DOS measured via Eq. ([5]) for tp > ipj. Typ- 
ical data are shown in Fig. [7] together with T = result 
obtained from the diagonalization of dynamical matrix in 
the ground state at T = 0. We observe that by increas- 
ing the temperature, deviations from the T = result 
emerge and become more pronounced. Notice, in partic- 
ular, that it is near the high-frequency cutoff that devi- 
ations appear first when T increases. By normalization, 
the low-frequency part of the DOS is then also affected. 
Qualitatively similar results were shown in Ref. 

For the particular density shown in Fig. [7J we observe 
that the DOS converges to its T — > limit when T < 
T* ps 1(T 7 . Thus, by contrast to the results of Ref. H3 
we conclude that a finite temperature region exists where 
the harmonic approximation holds. 

To quantify the size of this harmonic regime, we mea- 
sure the changes in d(u) through its first moment: 



u>i(tp,T) 



dio d(ui)u>. 



(30) 



The observed change in d(w) suggest that ui\ should de- 
crease with increasing T, because the DOS gives less 
weight to large frequencies when T increases. 

Our numerical analysis fully confirms this expectation, 
see Fig. [7J We find that, for any given density ip > ipj, 
the first moment wi decreases above some crossover tem- 
perature T* , which depends strongly on the volume frac- 
tion. To see this, we rescale in the lower panel of Fig. [7J 
the first moment lj\ by its value at T = 0, and present 
its evolution as a function of the rescaled temperature 
T/T*. Our numerical results are consistent with a sim- 
ple behavior 
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FIG. 7: Limit of validity of harmonic approximation above 
jamming. Top: The density of state d(uj) at ip = 0.655 
(<p—ipj ps 0.0084) measured for various temperatures through 
Eq. ([5}. The T = DOS is obtained from diagonalization of 
the dynamical matrix. Bottom: Convergence of the first mo- 
ment of the density of state to its T = limit below a crossover 
temperature scale T* = 10~ 3 ((p — <pj) 2 which is independent 
of system size. The vertical line indicates T/T* — 1. 



where W(x — > 0) = 1, with 

T* ps a (<p - <pjf, 



(32) 



wi(p,T = 0) \T* 



(31) 



where a ~ 10 -3 is a numerical prefactor. There is an am- 
biguity in the absolute value of the numerical coefficient 
a, since the scaling collapse in Fig. [7J does not depend on 
it. We have determined a such that W(x) starts to show 
deviations from 1 when x ss 1. 

To fully demonstrate that a harmonic regime exists 
in the thermodynamic limit, we must make sure that 
T* does not decrease to zero as the system size N is 
increased. Indeed, a strong system size dependence of 
the onset of anharmonicity was reported in Ref. 22( . To 
test this idea, we have repeated the above measurements 
for various system sizes, N = 120, 1000, and 8000. As 
shown in Fig.[7Jwc find that our results are devoid of any 
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FIG. 8: The structure of the three scaling regimes where 
power law divergences and 'anomalous' vibrational motion is 
observed. Regimes I and II are described by T = harmonic 
theories, while dynamics in Regime III is fully anharmonic. 
These regimes are separated by the crossover temperature 
T* ~ KT 3 ^ - ipj) 2 discussed in Sec. [VTU Note the small 
range of parameters where the effect of the T = jamming 
critical point are felt and Eqs. (|34l I37[) hold. 



system size dependence, although the statistical noise is 
of course larger for smaller systems. 

The fact that T* vanishes quadratically with (92 — ipj) 
for all system sizes confirms that materials exactly at 
ip j are 'marginally' solid with special properties. In par- 
ticular, the linear response and harmonic approximation 
regimes have a vanishing domain of validity [f| |45| . Our 
rcsults seem consistent with the analysis of low-energy 
barriers performed numerically in Ref. [||. In particular, 
we conclude that the T = DOS cannot be used at ipj 
to infer properties of the vibrational dynamics at any fi- 
nite temperature. This resolves the problem raised in the 
introduction. The material does not 'melt' at jamming 
when thermal fluctuations are added, because a finite 
amount of thermal noise pushes the system in a anhar- 
monic regime preventing the system from breaking apart. 

In conclusion, we find that for ip > ipj and T < T*(ip) 
the density of states converges to its zero temperature 
limit which is correctly described by the harmonic ap- 
proximation. In this regime the critical behavior is the 
one discussed previously in the zero temperature limit. 
We shall call this part of the phase diagram 'Regime I', 
see Fig|8l Above the temperature scale T*(ip) the har- 
monic approximation breaks down and the critical be- 
havior crosses over to Regime III that we discuss later. 
An explanation for the quadratic scaling of T*(tp) is pre- 
sented in Sec. IVII1 



B. Regime II: Effective harmonic description below 
jamming and its domain of validity 

Below jamming, the T — > limit inherently produces 
a strongly anharmonic system, because it produces hard 
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FIG. 9: Limit of validity of (effective) harmonic approxima- 
tion below jamming. Top: Temperature-scaled density of 
state d(uj)y/T at tp = 0.630 and various temperatures con- 
verges to T = limit below T* « 10 -7 . Bottom: Convergence 
of the height of the first peak of the DOS to its T = value be- 
low the crossover temperature scale T*(ip) = 10~ 3 (ipj — ip) 2 . 
The vertical line indicates T/T* — 1. 



sphere packings with zero energy, and the Hessian is not 
defined. However, an effective harmonic description that 
describes this behavior for tp < tpj can be approximately 
obtained Q. In this section we repeat the analysis per- 
formed in Regime I to determine the limit of validity of 
the effective harmonic description below jamming, which 
defines 'Regime IF in the phase diagram, see Fig. [8] 

In Fig. |9l we characterize the temperature evolution of 
the DOS at fixed volume fraction and decreasing temper- 
ature at ip = 0.63 < ipj. To remove a trivial temperature 
dependence of the time scale for hard spheres, we plotted 
d(uj)/VT versus luVT, as it is this rescaled DOS which 
eventually becomes independent of T as T — >• 0. 

Again, we observe that the DOS actually converges 
to its T — > limit for this particular density below 
T* rj 10~ 7 . To quantify the extent of the hard sphere 
regime, we focus on the height of the first peak of the 
DOS, di(T), which is a very sensitive indicator of the 
temperature evolution of d(ui). As in the case of the first 



13 



moment above jamming, we follow the convergence d\ (T) 
to its T = value, and determine the crossover tempera- 
ture T*(<p) controlling this convergence, see Fig. [9l The 
numerical results show that there is a again a simple be- 
havior for this quantity, 



in order to recover the correct T 



divergences 



di(<p,T = 0)~ \T* 



(33) 



where T>(x 0) = 1 with T* taken from the previous 
section. As in the case of jammed phase, there is a slight 
ambiguity in the actual prefactor of T*. For the un- 
jammed region, we simply use the a value determined for 
the jammed region, and find that this works very well. 

In conclusion, we find that for tp < tpj and T < T*(tp) 
the density of states converges to its zero temperature 
limit which is correctly described by the effective har- 
monic approximation [8j. In this Regime II, the criti- 
cal behavior is the one discussed previously in the zero 
temperature limit, see Fig[T2J For T > T*((p) the har- 
monic approximation breaks down and the critical be- 
havior crosses over to Regime III that will be discussed 
in the following. An explanation for the quadratic scaling 
of T*(<p) will be presented in Sec. IVHI 



C. Regime III: Anharmonic critical regime at finite 
temperature 

We found in previous sections that below the tempera- 
ture scale T*(tp), which vanishes quadratically and sym- 
metrically in \tp— ipj \ approaching the jamming transition 
from both sides, the critical behavior is well described by 
the harmonic theory, which leads to symmetric diver- 
gences on both sides of tpj, captured by 



X4 



\<p - v.j 



-1/2 



(34) 



We also showed that a finite temperature cuts off these 
divergences. This means that at tp = tpj all these observ- 
ables arc finite but diverge as T — ¥ 0. This represents 
a distinct critical regime which cannot be described by 
the harmonic approximation, and which is therefore con- 
trolled by anharmonicity and finite temperatures. We 
call this 'Regime IIP, see Fig [5] 

In order to obtain scaling laws in this regime, we use 
scaling arguments and assume that the behavior of the 
physical observables on both sides of the jamming tran- 
sition is connected via (observable-dependent) scaling 
functions 01 HI: 



0(<p,T) 



1 



\tp-tpj\V* 



H 



<P-VJ 



(35) 



where O stands for any of the following observables: Aoo, 
t, X4 or £|. The three regimes I, II, and III described 
in the phase diagram in Fig [12] correspond respectively 
to a; > 1, i < -1, and \x\ <C 1. Therefore, the scaling 
functions must be such that: 



Ho(x) 



Oi 



±00, 



(36) 



characteristic of regimes I and II, respectively; are 
observable-dependent numerical constants. 

The critical behavior in regime III can then be ob- 
tained by requiring that when tp — > tpj at small but finite 
temperature the dependence on tp — ipj drops out from 
Eq. (|35[) . In this way, one finds that at tp = tpj: 



Ac 



X4 



8~T 



-1/4 



(37) 



The critical behaviors of the three regimes are summa- 
rized in Fig. [8] Note that contrary to regimes I and II, 
physical or microscopic interpretation for the divergences 
in Eq. ([3"T]) are missing. A possible path could be the ex- 
tension of the effective harmonic description of Ref. Q 
to continous potentials, but this is beyond the scope of 
this work. 



VII. ANHARMONICITY AND THE 
TEMPERATURE SCALE T* 

We have found by numerical simulations that the har- 
monic description for the jammed and unjammed phase 
breaks down below a characteristic temperature scale 
vanishing quadratically with the distance to tpj. In the 
following we explore the origin of the anharmonicity re- 
sponsible for this behavior and offer an explanation for 
the scaling of T* with \tp — tpj\. 



Perturbative analysis of high-order 
nonlinearities 



A first natural idea to explain the breakdown of the 
harmonic approximation is to realize that the harmonic 
analysis corresponds to a truncation to second order of 
an expansion of the energy function around a global min- 
imum, which is generally justified if the temperature is 
low enough. Thus, the emergence of anharmonicity at 
finite temperature could be due to higher order terms in 
the expansion becoming relevant. 

To explore this hypothesis, we focus on the jammed 
phase where a genuine harmonic approximation holds, 
and an expansion of the energy can be performed. To 
this end, we compute the first low-temperature correc- 
tions to the average energy due to higher-order terms, 
and determine the temperature above which these non- 
linearities can no longer be neglected. 

We expand the position vector of a particle using a 
normal mode decomposition defined from a given energy 
minimum. We write 



(38) 



where rp' is the position of particle in the ground state 
and n a .i is the displacement of the z-th particle in the a- 
th normal mode. Using these normal modes, we can also 
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FIG. 10: The crossover temperature scale for the onset of an- 
harmonicity T* (bottom data with filled circles) compared to 
the temperature scale T pt (top data with open symbols differ- 
ent symbols for different packings) quantifying the relevance 
of the high order terms in the perturvative expansion of the 

1/2 

energy. The lines indicate T oc _E gs and T pt oc E g i . Clearly, 
these two temperature scales differ. 



sizes, but we have numerically checked that system size 
dependence is negligible, or at least much too weak to in- 
fluence our conclusions. For each configuration, we deter- 
mine numerically the normal modes at T = for various 
densities above jamming. From the numerical determi- 
nation of the modes, we estimate E\ and E2 in Eqs. (|40l 
|4"T} . and we deduce T pt from Eq. (|4"2"T) . We present our 
numerical results for T pt in Fig. [TU] where it is represented 
paramctrically as a function of the ground state energy, 
E gs (x [ip — ipj) 2 . For comparison, we also represent the 
crossover scale T* directly deduced from the study of the 
density of states d(uj), as discussed above. 

The numerical results indicate that T p t is larger than 
T* at all investigated volume fractions by several orders 
of magnitude, typically more than a factor 10 . Also, 
the density dependence of T pt is quite different from 
the one of T*, since the numerical data seem to indi- 
cate that T pt ~ (Egs/N) 1 / 2 ~ (tp — ipj). This means 
that the harmonic-anharmonic transition observed in the 
DOS cannot be explained through a pcrturbative anal- 
ysis around the ground state. These results show that 
thermally excited packings near jamming are not stabi- 
lized by quartic terms in the expansion of the energy, 
and another explanation has therefore to be sought for 
the existence of T*. 



expand the potential energy around the ground state en- 
ergy into a power series of {x a }. Taking the thermal 
average of this series, we obtain the following tempera- 
ture expansion form of the potential energy (details of 
this calculation are reported in the Appendix): 



(E) = E gs + E X T + E 2 T 2 + C(T 3 ) 
where E gs is the ground state energy, and 
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From these results, we can estimate a characteristic tem- 
perature, T p t, at which the quadratic term in the tem- 
perature expansion becomes equal to the linear one, sig- 
nalling that a harmonic expansion would break down: 



T pt = 



Ei 

En 



(42) 



We have computed T pt for 10 independent configura- 
tions with N = 256 particles. Due to the multiple sums 
in Eq. (|4ip we are forced to use relatively small system 



B. 'Fragility' of contact network and 
non-analyticity of pair potential 

The above calculation is pcrturbative in nature. It 
amounts to assuming that when heated, the system will 
simply explore the energy minimum used to compute the 
dynamical matrix. Crucially, this reasoning implicitly 
assumes that the forces are analytic functions, which is 
not correct for the harmonic pair potential we use which 
is not analytic in r = a where it is truncated. In fact, 
the above calculation would hold if v(r) = (1 — r/a) 2 
for all r- values. Indeed, it was found in Ref. [22| that the 
threshold for anharmonicity for the non-truncated poten- 
tial is about 10 4 larger than for the truncated potential. 
This numerical finding is in fact in quantitative agree- 
ment with the results of the previous section, recall the 
data shown in Fig. [TU1 

Therefore, we follow Ref. [22j and explore the idea 
that the change in the DOS observed above T* is as- 
sociated with the breaking and renewal of the 'contact 
network'. In such case, the anharmonicity would di- 
rectly result from the non-analyticity in the pair poten- 
tial. More precisely, it could be that some contacts are 
lost (or created) when T increases long before the high- 
order non-linearities become relevant. This physics can- 
not be described within a pcrturbative expansion around 
the ground state and is thus nonperturbativc in nature. 

In order to get a quantitative estimate of T* we must 
determine the temperature scale for which thermal fluc- 
tuations become large enough to disrupt the contact net- 
work. We start again with the jammed phase. Our phys- 
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ical idea is that since overlaps between particles are small 
close to jamming, they are more rapidly blurred by the 
effect of thermal vibrations. 

Let us now make this idea more precise. At volume 
fraction ip > ipj and T = 0, the typical length of the 
overlap between contacting particles is So = a— |rj — r^j ~ 
(ip — ipj). The latter estimate comes from the scaling of 
the width of the first peak in the pair correlation func- 
tion near ipj. On the other hand, at finite temperature, 
the overlap between two particles becomes a fluctuating 
quantity. The typical displacement between two neigh- 
boring particles has different scaling depending on the 
direction, longitudinal (L) or transverse (T), to the bond 



vector r» 



Sl oc Vt, St oc \/T(tp — ipj) 



-1/4 



(43) 



In order to loose a contact one must have either Sl ~ Sq 
or St ~ V~So, from purely geometrical reasoning. The 
first identity is more constraining, and implies that the 
notion of 'contacts' between particles only makes sense 
when y/T < (ip — ipj), which directly explains the scaling 
of the crossover temperature T* ~ (ip — </?j) 2 . 

On the hard sphere side below jamming, a similar rea- 
soning holds with the overlap being replaced by the gap 
between the particles. Note that this argument also im- 
plies that the quadratic dependence of T* with the dis- 
tance to if j is actually a direct consequence of the expo- 
nent in the pairwise interaction potential, which we have 
considered to be purely harmonic throughout this study. 

We can get an even more quantitative grasp of the 
crossover scale T* through an analysis of the radial dis- 
tribution function g(r). At T = 0, g(r) has a sharp peak 
near r = a, and its width is So — 0A(ip — ipj), the prefac- 
tor being determined numerically. In general, the radial 
distribution function is a volume average of the radial 
distribution function around a specific particle, 



fl(r) = ^E«( r ) = ^E(E*( p -i'«i) 



(44) 



At T = 0, gi(r) is the sum of delta functions correspond- 
ing to the relative distance between particle i and all its 
neighbors in the studied packing. When T > 0, these 
delta functions broaden under the influence of thermal 
fluctuations, and thus they acquire a width of the order 
s/T discussed above. 

In Fig. 111! we presents numerical results for g(r) and 
the gi(r) for a randomly chosen particle i at various tem- 
peratures both above and below T*, which for this ex- 
ample is T* w 4 • 10" 7 . At T = 10~ 8 , &(r) shows 
well-resolved delta functions, while g(r) is close to its 
T = value. With increasing temperature, the peaks 
in <7i(r) broaden, and cannot be resolved anymore above 
T* when their width becomes comparable to the width 
of the T = width of g(r). At even higher temperature, 
T = 10 -5 , the overall shape of g(r) is now controlled by 
the broadening of individual gi(r), and its shape becomes 
very different from the structure at T = 0, see Ref. [20j | . 
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FIG. 11: The radial distribution function g(r) at finite T, 
(black / thick lines), at T = (blue / thin line), and its 
particle-resolved version <?i(r) (red / dashed lines) for a ran- 
domly chosen particle i at ip = 0.66 and various temperatures. 
While contacts between i and its neighbors are well-defined at 
low T, they cannot be resolved when T becomes larger than 
T* which is w 4 x 10 -7 for this particular density. 



Let us push the argument further. Quantitatively, we 
determine numerically that the thermal broadening of in- 
dividual peaks in gt(r) is given by ~ 0.9y/T, such that 
at T = T* we obtain a width ~ 0.03(^ - ipj). This 
should be compared to the width of g(r) at half maxi- 
mum, So ~ 0A(ip — ipj), which is about 10 times larger. 
This numerical factor is reasonable since each particle 
possesses on average about 6 neighbors, and so the small- 
est of these overlaps is about Sq/6 ~ 0.066(</2 — ipj) which 
is indeed not far from 0.03(<y9 — ipj). 

Our argument suggests that anharmonicy starts to 
play a role and modify the DOS when an extensive num- 
ber of contacts become affected by thermal fluctuations. 
This is at odds with the results in Ref. [22j where a 
threshold for anharmonicity was defined via the breaking 
of a single contact in a given packing. This qualitative 
difference in our analysis presumably explains the differ- 
ent conclusion that we reach in this work, namely that 
a finite regime exists where a harmonic approximation is 
valid even in the thermodynamic limit, thus contradict- 
ing the different claim made in Ref. [22( that the domain 
of validity of harmonic theory is null. 

To summarize, we find that the onset of anharmonic- 
ity does not follow from the breakdown of a perturbative 
expansion to second order, but rather stems from the 
non-analycity of the pair interaction. Remark that the 
use of a truncated pair potential is of course a mandatory 
ingredient for the system to display a jamming transition 
in the first place. There is thus no room, it would seem, 
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to make the physics associated to the jamming transition 
more robust against thermal fluctuations since the very 
origin of the critical fluctuations — a marginally stable 
solid with vanishingly small overlaps between particles — 
also makes it highly sensitive to thermal fluctuations. 
This conclusion suggests that experimental investigations 
of the jamming transition using colloidal particles are not 
straightforward, as we now discuss. 



VIII. DISCUSSION: EXPERIMENTS VERSUS 
THEORY 

The theoretical analysis conducted in this article shows 
that single particle dynamics is able to directly reveal the 
divergence of time scales and length scales when the jam- 
ming transition at (T = 0, ip = ipj) is approached. We 
have additionnally shown that dynamic criticality sur- 
vives a finite amount of temperature, and have carefully 
explained under what conditions the signatures of the 
jamming singularity remain observable in the dynamics 
at finite temperature, emphasizing in particular the ex- 
istence of three distinct critical regimes in the vicinity of 
the critical point, Fig. [8] We suggested a possible connec- 
tions between the vibrational heterogeneity unveiled here 



and the experimental reports in Rcf. [26J, |27j . Other ex- 
periments report a nonmonotonic evolution of dynamic 
heterogeneity with density [48|, [4{|, but the connection 
with our work is less clear. 

We mentioned in the introduction the important ex- 
perimental activity in the colloidal community aiming 
at characterizing the emergence of 'soft modes' near the 
jamming transition. It is interesting to use our findings 
to discuss these experimental investigations. Two differ- 
ent types of colloidal systems have been used, which we 
discuss separately. 

A first series of experiments concentrates on colloidal 
PMMA hard spheres @, H ■ In that case, temper- 
ature simply serves to produce Brownian motion but 
since particles are 'infinitely' hard, temperature docs 
not blur the jamming transition, and thus these exper- 
iments arc in principle well-suited to probe Regime II, 
the critical regime of hard spheres. However, the inves- 
tigations in Refs. [H [l(| were performed in the range 
p = 0.56 — 0.60. Although some care must be taken 
with the absolute values of p> in experimental work, 
these values seem consistent because the amplitude of 
the experimentally-measured DW factors are in the range 
A 2 (oo) ps (0.005 - 0.03)cr 2 , see Fig. [TJ However, these 
volume fractions are very far from pj (see Figs. 2] and 
[5]) and all lie outside the critical regime II where a sharp 
drop in the DW factor would be observed, accompanied 
by the appearance of low-frequency modes and a grow- 
ing correlation length. It would therefore be interest- 
ing to push these experiments further towards the jam- 
ming transition, to see whether experimental signatures 
of the dynamic criticality discussed in this work can be 
detected. In Fig. 1121 we provide an extended view of 
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FIG. 12: An enlarged view of the (Temperature, Volume frac- 
tion) phase diagram reporting the critical regimes shown in 
Fig. [SI and the approximate location of the experimental stud- 
ies aimed at studying anomalous vibrational dynamics in col- 
loidal systems: Gosh et al. [i| use PMMA hard spheres, while 
Chen et al. [IH and Caswell et al. [l5[ study PNIPAM micro- 
gel particles. Previous studies all lie too far away from the 
jamming transition to detect the dynamic criticality associ- 
ated to the transition. 



the (T, ip) phase diagram including the experimentally 
explored regime for hard spheres, outside the critical re- 
gion discussed in our work. 

A second series of experiments designed to probe low- 
frequency modes near jamming uses soft 'PNIPAM' mi- 
crogel particles flll - [l3[ [Pol ] . In this case both the vol- 
ume fraction and the temperature, expressed in units 
of the particle softness, are relevant control parameters. 
Because microgel particles arc easily compressed, the 
volume fraction can be adjusted to be as close to the 
jamming density as desired. However, softness is now 
the main issue. The adimensional temperature estimate 
is T/e ps 2 x 10~ 5 in Ref. [Tl|. We can confirm this 
value by again using the DW factor as a sensitive 'ther- 
mometer'. Indeed, the values measured near jamming, 
A 2 (oo) ~ 3 x 10 _4 (T 2 , correspond in our Fig. [J to a 
temperature scale of about 10 -5 . In very recent work 
with different microgel particles [lj|, the DW factor is 
larger, A 2 (oo) ~ 2 x 10 -3 cr 2 , and thus the system is 
effectively at an even higher temperature, of the order 
10 -4 . On the basis of our numerical results, we conclude 
that the microgel particles studied in Refs. [IJ EH are 
so soft (or, cquivalcntly so 'hot') that anharmonicity is 
very strong because temperature is much larger than the 
crossover scale T*, see their location in the phase dia- 
gram of Fig. [12] As for hard colloids, we suggest that it 
could be interesting to design different experimental sys- 
tems with particles that are less soft than microgels, for 
instance emulsions, in order to probe experimentally the 
critical dynamics associated to the jamming transition. 

In conclusion, because of the smallness of the critical 
region (that certainly would deserve a theoretical expla- 
nation) experiments have not yet probed critical prop- 
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erties related to the jamming transition. Since they are 
quite far from the critical region, our results imply that 
there cannot be any direct connection between the DOS 
measured in those experiments and the theoretical re- 
sults established at T = 0. On the basis of our results, 
one should find an alternative explanation, not based on 
'thermal vestiges' of the jamming transition, for the ex- 
perimental results. Indeed, there is no need to invoke any 
critical behavior to account for the smooth evolution with 
density of the low- frequency uj* reported in Ref. [ll[ . It is 
fully consistent with the data obtained in Fig. [5] for large 
temperatures which basically track the microscopic time 
scale To that changes smoothly across ipj as a direct re- 
sult of the compression. Another experimental signature 
of the jamming singularity is the nonmonotonic evolution 
of the pair correlation upon compression [l5l . l47j . This 
can also be explained in a way unrelated to the critical 
properties of t he j amming transition since, as demon- 
strated in Ref. [l9(, this anomaly persists arbitrarily far 
from the transition, even in the equilibrium fluid. 

In conclusion, our study reveals that single particle mo- 
tion in thermal systems near jamming becomes singular, 
and reveals a number of scaling laws and divergences for 
both time scales and length scales. This dynamic criti- 
cality thus provides direct evidence of the critical nature 
of the jamming transition occurring at T = 0, and might 
suggest new routes to attack the problem from a theo- 
retical viewpoint, for instance using the replica approach 
in the framework of statistical mechanics. In addition, 
we have suggested that it would be very interesting to 
design new experiments using colloidal systems to probe 
more closely the diverging time scales and length scales 
associated to vibrational dynamics of dense packings. 
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Appendix A: Derivation of Eqs. (|39l- I41|) 

Using the normal modes decomposition of the position 
vector of a particle Eq. (|38|). the potential energy of the 
system can be expanded as 

E = E gs + i ^ VabXaXb + ^ ^ VabcXaXbXc 
ab abc 

+ 0(x 5 ), (Al) 

abed 

where we introduced a simplified notation for the deriva- 
tives, e.g. V a b = g° g Xb ■ Note that the first order term 
vanishes since we focus on a potential energy minimum. 
We regard this development of the potential energy as 
the sum of the harmonic term Eq (up to second order) 
and anharmonic terms AE. When AE is treated as per- 
turbation, we can obtain the following expansion of the 
distribution function, 

/ = f (l-f3(AE-(AE)o) + l/2f3 2 (AE 2 -(AE 2 ) Q 
-2AE(AE) + 2{AE) 2 ) + G(/3 3 A£ 3 )) , (A2) 

where / is the canonical distribution function corre- 
sponding to E and /□ to Eq; (■ ■ -)o indicates a thermal 
average over the distribution /q. Taking the thermal 
average of Eq. (|Al|) using the distribution function in 
Eq. (|A2[) . we obtain the perturbation expansion of the 
potential energy. 

To organize this series as the low-temperature expan- 
sion, we have estimated the temperature dependence of 
each term using the fact that (x n )o cx T n l 2 when n is 
even. As a result, we get the following equation: 



(E) = E gs + - ^ Vab(x a Xb)o + X! V a.bcd(x a XbX c Xd)o - ^ (jjjgj V a bcVdef(x a X b X c X d X e Xf} 

ab abed abcdef 

+ 7^V ab V c def(x a X b X c X d X e X f ) ~ -^y ab V c d e f{x a X b ) (x c XdX d X f ) )^ 

+ 2^2 X! -^^V ab V cde Vfgh({XaXbX c XdXeXfXgXh)o- (x a X b )o(x c XdX e XfXgXh)o^ +0(T 3 ), (A3) 
abcdef gh 

I 

where the first term is constant, the second term is pro- 
portional to T, the next 6 terms are 0(T 2 ), and the rest 
is 0(T 3 ). Although four-, six- and eight-point functions 
appear in this expression, they can be evaluated using 
the Gaussian approximation and Wick's theorem. For 



example, the expectation value of the two-point function 
is simply 

(x a x b )o = (A4) 
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where A a is the eigenvalue for the a-th normal modes. with Wick's theorem, we obtain the final expressions in 
Note also that the second order derivatives are the Eqs. p9H41[) . 
diagonal matrix V a t = ^ a o~ab- Using these relations 
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